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Abstract 

One- and two-parameter families of flows in H 3 near an Andronov-Hopf bifurcation (AHB) 
are investigated in this work. We identify conditions on the global vector field, which yield a 
rich family of multimodal orbits passing close to a weakly unstable saddle-focus and perform 
a detailed asymptotic analysis of the trajectories in the vicinity of the saddle-focus. Our anal- 
ysis covers both cases of sub- and supercritical AHB. For the supercritical case, we find that 
the periodic orbits born from the AHB are bimodal when viewed in the frame of coordinates 
generated by the linearization about the bifurcating equilibrium. If the AHB is subcritical, it 
is accompanied by the appearance of multimodal orbits, which consist of long series of nearly 
harmonic oscillations separated by large amplitude spikes. We analyze the dependence of the 
interspike intervals (which can be extremely long) on the control parameters. In particular, 
we show that the interspike intervals grow logarithmically as the boundary between regions of 
sub- and supercritical AHB is approached in the parameter space. We also identify a window 
of complex and possibly chaotic oscillations near the boundary between the regions of sub- 
and supercritical AHB and explain the mechanism generating these oscillations. This work is 
motivated by the numerical results for a finite-dimensional approximation of a free boundary 
problem modeling solid fuel combustion. 



1 Introduction 



Systems of differential equations in both finite- and infinite-dimensional settings close to an AHB 
have been subject to intense research due to their dynamical complexity and importance in ap- 
plications. The latter range from models in fluid dynamics [29J to those in the life sciences, in 
particular, in computational neuroscience [3], [TT1 [1^1 Ell [23 EQl EHl [36] . When the proximity to the 
AHB coincides with certain global properties of the vector field, it may result in a very complex 
dynamics [U [9j [TO] [2"T| [23] . The formation of the Smale horseshoes in systems with a homoclinic 
orbit to a saddle-focus equilibrium provides one of the most representative examples of this type 
[2]. Canard explosion in relaxation systems affords another example (6[ [27]. Recent studies of 
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relaxation systems, motivated mainly by applications in the life sciences, have revealed that the 
proximity to an AHB has a significant impact on the system dynamics. It manifests itself as a 
family of multimodal periodic solution that are composed of large-amplitude relaxation oscillations 
(generated by the global structure of the vector field) and small-amplitude nearly harmonic oscil- 
lations (generated by the vector field near the equilibrium close to the AHB) [T21 [Ml EH ED EH EZ] 
(see Figure [1]). These families of solutions possess rich bifurcation structure. A remarkable example 
of an infinite-dimensional system close to the AHB has been recently studied by Frankel and Royt- 
burd [TH [15j HU [T71 [18] . They derived and systematically studied a model of solid fuel combustion 
in the form of a free boundary problem for a ID heat equation with nonlinear conditions imposed 
at the free boundary modeling the interface between solid fuel mixture and a solid product. The 
investigations of this model revealed a wealth of spatial-temporal patterns ranging from a uniform 
front propagation to periodic and aperiodic front oscillations. The transitions between different 
dynamical regimes involve a variety of nontrivial bifurcation phenomena including period-doubling 
cascades, period-adding sequences, and windows of chaotic dynamics. To elucidate the mechanisms 
responsible for different dynamical regimes and transitions between them, Frankel and Roytburd 
employed pseudo-spectral techniques to derive a finite-dimensional approximation for the interface 
dynamics in the free boundary problem |17j . As shown in [13} 117] . a system of three ordinary 
differential equations captured the essential features of the bifurcation structure of the infinite- 
dimensional problem. The numerical bifurcation analysis of the finite-dimensional approximation 
revealed a rich family of multimodal periodic solutions similar to those reported in the context of 
relaxation systems near the AHB [13J. The bifurcation diagrams presented in |13j and in [30] share 
a striking similarity, despite the absence of any apparent common structures in the underlying 
models (except to the proximity to the AHB). In particular, in both models, topologically distinct 
multimodal periodic solutions are located on isolas, closed curves in the parameter space. The 
methods of analysis of the mixed-mode solutions in [261 1301 l3"T| [38] used in an essential way the 
relaxation structure present in these problems. These approaches can not be applied directly to 
analyzing the model in [13] , because it is not a priori clear what creates the separation of the time 
scales in this model, in spite of the evident fast-slow character of the numerical solutions. This is 
partly due to the spectral method, which was used to derive the system of equations in [13] : while 
it has captured well the finite-dimensional attr actor of the interface dynamics, it has disguised the 
structure of the physical model. One of the goals of the present paper is to identify the structure 
responsible for the generation of the multimodal oscillations in a finite-dimensional model for the 
interface dynamics and to relate it to those studied in the context of relaxation oscillations. 

The family of flows in [13] includes in a natural way two types of the AHBs. Depending on the 
parameter values, the equilibrium of the system of ordinary differential equations in [13] undergoes 
either a sub- or a supercritical AHB. A similar situation is encountered in certain neuronal models 
(see, e.g., [HI [33]). In either case, the global multimodal periodic solutions are created after the 
AHB. However, in the case of a supercritical bifurcation, they are preceded by a series of period- 
doubling bifurcations of small amplitude limit cycles, arising from the AHB. On the other hand, in 
the subcritical case, the AHB gives rise to multimodal solutions, whose lengths and time intervals 
between successive large amplitude oscillations can be very long. In the present paper, we perform 
a detailed asymptotic analysis of the trajectories in a class of systems motivated by the problem 
in [13]. Our analysis includes both cases of the sub- and supercritical AHBs. We also investigate 
the dynamical regimes arising near the border between the regions of sub- and supercritical AHB. 
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This region in the parameter space contains a number of nontrivial oscillatory patterns including 
multimodal trajectories with substantial time intervals between successive spikes, irregular, and 
possibly chaotic oscillations, as well as a family of periodic orbits undergoing a cascade of period- 
doubling bifurcations. Our analysis shows that these dynamical patterns and the order in which 
they appear under the variation of the control parameters are independent on the details of the 
model, but are characteristic to the transition from sub- to supercritical AHB. 

The outline of the paper is as follows. After introducing the model and rewriting it in the normal 
coordinates, we present a set of the numerical experiments to be explained in the remainder of the 
paper. Then we state our results for each of the following cases: supercritical AHB, subcritical 
AHB, and the transition layer between the regions of sub- and supercritical AHB. In Section 
3, we analyze the local behavior of trajectories near a weakly unstable saddle-focus. The local 
expansions used in this section are similar to those used in [8] for analyzing the AHB using the 
method of averaging. However, rather than establishing existence of periodic solutions, the goal of 
the present section is to approximate the trajectories near the saddle- focus. For this, after rescaling 
the variables and recasting the system into cylindrical coordinates, we reduce the system dynamics 
to a 2D slow manifold. By integrating the leading order approximation of the reduced system, we 
obtain necessary information about the local behavior of trajectories. The results of this section are 
summarized in Theorem 3.1. In Section 4, we study oscillatory patterns generated by the class of 
systems under investigation. It is divided into three subsections devoted to the oscillations triggered 
by the supercritical AHB (§4.1), subcritical AHB (§4.2), and those found in the transition region 
between sub- and supercritical AHB (§4.3). In the supercritical case, we show that the oscillations 
just after the AHB are already bimodal. Generically, one needs to use two harmonics to describe the 
limit cycle born at the supercritical AHB in IR 3 . We give a geometric interpretation of this effect, 
which we call the frequency doubling, due to the fact that the second frequency is twice as large 
as that predicted by the Hopf Bifurcation Theorem [29] . We also compute the curvature and the 
torsion of the periodic orbit as a curve in IR 3 . The latter is useful for the geometric explanation of 
the frequency doubling. In Subsection 4.2, we study the subcritical case. We show that under two 
general assumptions on the global vector field, the presence of the return mechanism (Gl) and the 
strong contraction property (G2), the subcritical AHB results in sustained multimodal oscillations. 
Even though the oscillations may not be periodic (our assumptions do not warrant periodicity), the 
time intervals between consecutive spikes of the resultant motion comply to the uniform bounds 
given in Theorem 4.1. This subsection also contains the definitions of the multimodal oscillations 
and the precise formulation of the assumptions on the global vector field. The proof of Theorem 
4.1 is relegated to Section 5. Subsection 4.3 describes the transition between sub- and supercritical 
AHB. This transition contains a distinct bifurcation scenario, when a small positive real part of 
the complex conjugate pair of eigenvalues is positive and fixed and the first Lyapunov coefficient 
changes sign. The numerical simulations show that this results in the formation of the chaotic 
attractor followed by the reverse period-doubling cascade. To explain this bifurcation sequence, 
we derive a ID first return map. The asymptotic analysis in this subsection is complemented by 
numerical extension of the map to the region inaccessible by local asymptotic expansions. The 
first return map obtained by the combination of the analytic and numerical techniques reveals 
the principal traits of the transition from sub- to supercritical AHB. Finally, the discussion of the 
results of the present paper and their relation to the previous work is given in Section 6. 
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2 The model and numerical results 



In the present section, we formulate the model and present a set of numerical results, which moti- 
vated our study. The following system of three ODEs was derived in |13| . as a finite-dimensional 
approximation for the interface dynamics in a free boundary problem modeling solid fuel combus- 
tion: 

3(v 3 + v 2 - vi) - vk(vi) - 

Vl = ' ( } 

v 2 = v 3 -vi, (2.2) 
v 3 = Q( Vl - V3 ) - Q V2 + u ( Vl + l)k( Vl ) + 2vf. (2.3) 

Here, v\{t) approximates the velocity of the interface between the solid fuel and the burnt mate- 
rial. Functions V2(t) and v 3 (t) are the first two coefficients in the series expansion for the spatial 
temperature profile with respect to the basis of Chebyshev-Laguerre polynomials. Equations (12. ip - 
(|2.3p are obtained by projecting the original infinite-dimensional problem onto a finite-dimensional 
function space and using the method of collocations for determining the unknown coefficients. An 
important ingredient of the model, nonlinear kinetic function k(v) is given by 

fl _ V )P - (l- 

k(v) = { - } p+ \ } . (2.4) 

It reflects the dependence of the velocity of propagation on the temperature in the front. This 
relation is very complex and is not completely understood at present. To account for a range of 
possible kinetic mechanisms, the model includes two control parameters: v and p. Both the interface 
velocity and the temperature profile are calculated in the uniformly moving frame of reference. 
Therefore, system of equations (|2.1|) - (j2.3[) describes the deviations of the interface dynamics from 
that of a front traveling with constant speed. In particular, periodic and aperiodic oscillations 
generated by (I2,ip - (l2,3p correspond to complex spatio-temporal patterns in the infinite-dimensional 
model. The numerical results presented in [13] show a remarkable similarity between the complex 
patterns generated by the infinite- and finite-dimensional models and between the scenarios for 
transitions between different regimes in both models. For more information about the derivation 
of the free-boundary problem for solid fuel combustion and its finite-dimensional approximation, 
we refer the reader to |13j and bibliography therein. 

A simple inspection of (|2.ip - (|2.3p shows that it has an equilibrium at the origin O = (0, 0, 0) T 
for all values of v and p > 0. We linearize (|2.ip - (|2.3p about the equilibrium at the origin 

%-v -3 

v v v \ T 
A(v)v + . . . , where A(v) = | — 1 1 |, v = (vi,V2,v 3 ) . 

9-u -6 

Near vah = \, Jacobian matrix A(y) has a negative eigenvalue —1 and a pair of complex conjugate 
eigenvalues A = a(v) + ib(u) and A. The latter crosses the imaginary axis transversally at v = i>ah'- 

a {vah) = and a' {vah) ^ 0. 
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Figure 1: Stable periodic solutions of system of differential equations (|2.5p . Plots in a-d show 
solutions near a supercritical AHB (p = 2.2) and those in e-h correspond to the subcritical case 
(p = 3). The time series of x\ are presented in the left column and those of X2 are given in the 
right column. The time series for xs and x% are qualitatively similar. 
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Figure 2: a. A periodic trajectory of the original system in the regime close 
to a subcritical AHB. Cross-sections S + and S" are defined as follows: S + = 
{(xi, pcos 9, psia.9) : p = 0.1, x\ S [—0.02,0.02], #€ [0, 27r)}. E~ is orthogonal to x\— axis and 
is set at x\ = —0.65. The numerical results shown in b test the return mechanism and the strong 
contraction property for (|2.5|) (see (Gl) and (G2) in the text). For this, we covered S + with a 
uniform mesh S^. N . = {Cij} , i = 1,2,..., Ni, j = 1, 2, . . . , Nj. Using points £ij as initial con- 
ditions, we integrated (|2.5p until the corresponding trajectories hit the second cross-section, E~, 
transverse to the stable manifold. The top plots in b shows the image of Sat^at m E~ under the 
flow-defined map. To test the strong contraction property, we repeated this experiment by taking 
E~ progressively closer to the origin (the location of S~ is indicated in each plot). The top five 
images in b clearly indicate that the projection of the vector field onto S~ is strongly contracting. 
The bottom plot shows that the trajectories starting from S + reach a small neighborhood of the 
unstable equilibrium. 



Therefore, the equilibrium of (|2.ip - (l2,3j) undergoes an AHB. In the neighborhood of vah, ot = a(v) 
defines a smooth invertible function. We shall use a as a new control parameter. After a linear 
coordinate transformation, system of equations (|2.ip - (|2,3p has the following form: 





x + h(x,a), x = (xx, x 2 , x 3 y , h = (h 1 ,h 2 ,h 3 ) T . (2.5) 



Here, by b(a) := b (a x (a)) we denote the imaginary part of A as a function of the new control 
parameter. Nonlinear function h : 1R 3 x IR — > 1R 3 is a smooth function such that h (0, a) = and 
^5: = f° r vames °f a close to 0. More precisely, h stands for a family of functions parametrized 
by p (see (|2.4p ). To keep the notation simple, we omit the dependence of h on the second control 
parameter p. It turns out that in the range of parameters of interest, the AHB of the equilibrium 
for a = can be either subcritical, or supercritical, or degenerate depending on the value of p. This 
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Figure 3: The duration of the ISIs is plotted for: (a) 7 = 0.778 and (b) a = 0.045. 



gives rise to several qualitatively distinct oscillatory regimes generated by (|2.5p for small values of 
a > 0. The corresponding bifurcation scenarios for (|2.ip - (|2.3p for fixed values of p and varying v 
were studied numerically in [13]. Below, we reproduce some of these numerics for the system in new 
coordinates and supplement them with a set of new numerical experiments relevant to the analysis 
of this paper. After that we state our results for each of the following cases: supercritical AHB, 
subcritical AHB, and the transition layer between the regions of sub- and supercritical bifurcations. 

We start with discussing the supercitical case. It is well known that a supercritical AHB 
produces a stable periodic orbit in a small neighborhood of the bifurcating equilibrium. The period 
of the nascent orbit is approximately equal to 2-7T/3 -1 , where j3 = 6(0) is the imaginary part of the 
pair of complex conjugate eigenvalues of the matrix of the linearized system at the bifurcation. The 
numerical time series of x\ and 0:2,3 show that, while X2.3 oscillate with the frequency prescribed by 
the Andronov-Hopf Bifurcation Theorem |28[ [29] , the former oscillates with twice that frequency 
(see Figure [T^,b). The asymptotic analysis of Section 3 explains this counter- intuitive effect and 
shows that this is, in fact, a generic property of the periodic orbits born from a supercritical 
bifurcation in R n , n > 3. Note that the difference in frequencies of oscillations in x\ and £2,3 can 
not be understood from the topological normal form from the AHB [lj, because the latter does not 
contain the information about the geometry of the periodic orbit. In Section 4.1, we compute two 
geometric invariants of the periodic orbits as a curve in IR 3 : the curvature and the torsion. The 
latter shows that generically the orbit born from a supercritical AHB is not planar. This together 
with certain symmetry properties of the orbit explains the frequency doubling of the oscillations 
in x\. As was noted in [13], for increasing values of a > 0, the small periodic orbit born from the 
supercritical AHB undergoes a cascade of period-doubling bifurcations, the first of which is shown 
in Figure [Tb,d. Already at the first period-doubling bifurcation, the periodic orbit lies outside the 
region of validity of the local power series expansions, developed in the present paper. Therefore, 
our analysis does not explain the period-doubling bifurcations for increasing values of a > 0. In 
Section 4.3 we complement our analytical results with the numerical construction of the ID first 
return map. The latter explains the mechanism for period-doubling cascade and the window of 
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Figure 4: The time series of X2 are plotted for decreasing values of 7 > 0. In plots a and b, 
the number of small amplitude oscillations separating the spikes increases, while the qualitative 
form of the solutions remains the same. For smaller values of 7, the system exhibits intermittency: 
series of regular small amplitude oscillations are followed by those of irregular oscillations of the 
intermediate amplitude (c and d). 



complex dynamics reported in [13j . 

We next turn to the subcritical case. The dynamics resulting from the subcritical AHB depends 
on the properties of the vector field outside of a small neighborhood of the equilibrium. We 
conducted a series of numerical experiments to study the global properties of the vector field (see 
caption of Figure [2] for details). Based on these numerical observations, we identify two essential 
features of the vector field: (Gl) the return mechanism and (G2) the strong contraction property. 
Specifically, 



(Gl) For a suitably chosen cylindrical crossection, S + , placed sufficiently close to the origin, 
and another crossection, E~ , transverse to W s (0) (see Figure [2K), the flow-defined map 
Q : E + — > E~ is well-defined for a > and depends smoothly on the parameters of the 
system. 

(G2) There is a region EE, adjacent to E~ and containing a subset of W s (0) (see Figure EK), 
in which the projection of the vector field in direction transverse to W s {0) is sufficiently 
stronger than that in the tangential direction. 



In Section 4.2, conditions (Gl) and (G2) are made precise. Properties (Gl) and (G2) guarantee 
that for small values of a > 0, the trajectories leaving a small neighborhood of the saddle-focus, 
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Figure 5: The time series of xi (left column) and the corresponding trajectories of (|2.ip - (l2.3p (right 
column) plotted for negative values of 7 close to 0. The transition from the complex dynamics in a,b 
to regular approximately harmonic oscillations in g,h contains a reverse cascade of period-doubling 
bifurcations (of). 
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after some relatively short time enter II. In II, the trajectories approach W s {0) closely and then 
follow it to a sufficiently small neighborhood of the unstable equilibrium. For small a > 0, the 
trajectories starting from a sufficiently small neighborhood of the saddle-focus, remain in some 
larger (but still small) neighborhood of the origin for a long time. Eventually, they hit S + and 
the dynamics described above repeats. Therefore, under conditions (Gl) and (G2), the subcrit- 
ical AHB bifurcation results in a sustained motion consisting of long intervals of time that the 
trajectory spends near a weakly unstable saddle-focus and relatively brief excursions outside of 
a small neighborhood of the origin (see Figure [2h). In terms of the time series, the dynamical 
variables undergo series of small amplitude oscillations alternating with large spikes (Figure [T^-h) . 
For increasing values of a > 0, the number of small amplitude oscillations decreases and so do 
the time intervals between consecutive spikes, so-called interspike intervals (ISIs). The multimodal 
solutions arising near the subcritical AHB are not necessarily periodic. Nonetheless, the ISIs may 
be characterized in terms of the control parameters present in the system, regardless of whether 
the underlying trajectories are periodic or not. Specifically, in Section 4.2, we derive an asymptotic 
relation for the duration of the ISIs: 



where the first Lyapunov coefficient 7(p) reflects the distance of the system from the transition 
from sub- to supercritical AHB, with j(p) being positive when the AHB is subcritical and equal 
to zero at the transition point. Finally, e > is a small parameter and positive a = 0(e 2 ). 
The role of e will become clear later. We illustrate (12. 6p with numerically computed plots of the 
ISIs under the variation of control parameters in Figure [3j The ISIs provide a convenient and 
important characteristics of the oscillatory patterns involving pronounced spikes. For example, 
it is widely used in both theoretical and experimental neuroscience for description of patterns of 
electric activity in neural cells. An important aspect of the ISIs, is that the dependence of the ISIs 
on the control parameters in the system (which often can be directly established experimentally) 
reveals the bifurcation structure of the system. Therefore, the analytical characterization of the 
ISIs in terms of the bifurcation parameters, such as given in (|2.6p . is important in applications. For 
the problem at hand, the ISIs depend on the interplay of the two parameters a and j(p), which 
reflect the proximity of the system to a codimension 2 bifurcation. The dependence of the ISIs 
on a was studied in |21j for a model problem near a Hopf-homoclinic bifurcation. Our analysis 
extends the formula for the ISIs obtained in [21] to a wide class of systems and emphasizes the 
role of "f(p), which is important for the present problem. We also note that the structure of the 
global vector field suggests that (|2.5p is also close to a homoclinic bifurcation, which can also 
influence the ISIs. However, under the variation of the control parameters a and p the system 
remains bounded away from the homoclinic bifurcation, so that the latter effectively does not 
affect the ISIs. In fact, (j2.6|) can be easily extended to include the distance from the homoclinic 
bifurcation as well. As follows from (12. 6ft . the ISIs increase for decreasing values of 7 > (see 
Figure [3b). This observation prompted our interest in investigating the transition in the system 
dynamics as 7 crosses 0. Numerical simulations show that as 7 approaches from the positive 
side, the oscillations become very irregular and very likely to be chaotic (Figure H]). When a > 
is fixed and 7 crosses 0, the irregular dynamics is followed by the reverse period-doubling cascade 
and terminates with the creation of a regular small amplitude periodic orbit, which can be then 
tracked down to a nondegenerate supercritical AHB (Figure (5b.-h) . This scenario presents a certain 




(2.6) 
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interest as it suggests the formation of the chaotic attractor in a well-defined bifurcation setting. 
Namely, the transition from the region of subcritical to supercritical AHB for small a > leads 
to the formation of a chaotic attractor and a reverse period-doubling cascade. As in the case of 
the complex dynamics appearing for increasing values of a for supercritical AHB, the irregular 
oscillations in the present case can not be understood using the local analysis alone. Therefore, 
we complemented our analysis with the study of numerically constructed first return map. The 
latter gives a clear geometric picture of the origins of the chaotic dynamics and the period-doubling 
cascade of periodic orbits in this parameter regime. 



3 The local expansions 

In the present section, we study trajectories of (12. 5p in a small neighborhood of the equilibrium at 
the origin. Assume that h(x,a) is a smooth function in a small neighborhood of (0,0), so that all 
power series expansions below are justified. We are interested in the dynamics of (12. 5p for values 
of a close to 0. The role of a in our analysis is twofold. On the one hand, a is a control parameter 
in this problem, on the other hand, the smallness of a is used in the asymptotic analysis below. To 
separate these two roles, we use the following rescaling 

a = pe 2 , (3-1) 

where e > is a fixed sufficiently small constant and p varies in a certain interval of size 0(1). We 
introduce cylindrical coordinates in IR 3 

(xi,X2,X3) i ^ (xi,p,0) eEx ]R + xS 1 , X2 = pcosO, X3 = psinO, (3-2) 

and define 

D = {(xi,p, 9) : \x 1 \<M 2 e 2 ,pe(0,2Me]},D = {(x 1 ,p,9): \x±\ < M 2 e 2 , p G (0, Me] } , (3.3) 
where M > is sufficiently large. 

To characterize the trajectories of (|2.5p in D we will use an exponentially stable slow manifold, 
S, whose leading order approximation is given by 

S = {(x 1} p,9) : xi = U(9)p 2 , 0<P< 2eM} , (3.4) 

where 

U(9) = a + Acos(29-$). (3.5) 
and a, A and $ are computable constants (see Appendix A). 

Below we show that the trajectories of (|2.5[) with initial data from Dq approach an O(e) neigh- 
borhood of So in time O ( | In e | ) and stay in this neighborhood as long as they remain in D. The 
reduction of the system dynamics to S yields a complete description of the trajectories of (|2.5p in D. 
The qualitative character of the solution behavior in D depends on the sign of the first Lyapunov 
coefficient 

7 = ^ J* {U (9) Q 2 (9) + 3 (9)) d9, (3.6) 



11 



where 2tt— periodic trigonometric polynomial Q 2 ,3 are given in Appendix A. Finally, by (3 = 6(0) 
we denote the absolute value of the imaginary parts of the complex conjugate pair of eigenvalues 
at the AHB (see (|2.5p ). The results of the analysis of this section are summarized in 

Theorem 3.1 Suppose (3 > and 7 ^ 0. Let e > denote a small parameter. Then 

FOR a = 0(e 2 ) AND FOR SUFFICIENTLY LARGE (INDEPENDENT OF e) M > 0, THE TRAJECTORIES 
WITH INITIAL CONDITIONS IN D ENTER AN O(e) NEIGHBORHOOD OF So IN TIME 0(|lne|) AND 
REMAIN THERE AS LONG AS THEY STAY IN D. IN THE NEIGHBORHOOD OF So, THE TRAJECTORIES 
CAN BE UNIFORMLY APPROXIMATED ON ANY INTERVAL OF TIME [t ,t\ OF LENGTH O(l): 

Xl (t) = p 2 U(8) + 0(e 3 ), x 2 (t) = pcos9 + 0(e 2 ), and x 3 {t) = psm9 + 0(e 2 ), (3.7) 

WHERE U(9) AND 7 ARE DEFINED IN (13. 5 D AND (13.61) . RESPECTIVELY, AND 

P (t) = ((p(t )- 2 + l)e- 2 ^-l)^+O(e% 9 = (3, tan(0(O)) = ^. 
\\ a/ a/ xi in 



Remark 3.1 



a) From Theorem 3.1, one can easily deduce the existence of a periodic orbit O a , when 07 < 0. 

O a is stable (unstable) if a > (a < 0). The leading order approximation for O a follows 
from (pO) : 

x x = p 2 U{9) +C»(e 3 ), x 2 = pcos(9 + 0(e 2 ), x 3 = psin 6 + 0(e 2 ), and 9 € [0, 2vr), (3.8) 

where p = <jl^- Equations (13. 5p and f)3.8|) imply that the frequency of oscillations in x\ is 

twice that of oscillations in ^2,3, provided A ^ in (|3 .51) (see Figure d^,b) . In Section 4.1, we 
give a geometric interpretation of this frequency doubling effect. 

b) If a > and 7 > 0, Theorem 3.1 describes the trajectories with initial conditions near a weakly 

unstable equilibrium. It shows that along these trajectories, ^2,3 undergo approximately 
harmonic oscillations, whose amplitude grows at the rate 0(e 2 ). In addition, they satisfy the 
following scaling relation: 

xi (i) ~ ap 2 (t), (3.9) 
where xi(t) stands for the average value of x\(t) over one cycle of oscillations. 

c) The domain of validity of the asymptotic expansions in (|3.7p extends much farther than D. 

As will follow from the proof of Theorem 3.1, M in the definition of D may be taken up to 

o fe~3~J. This implies that, with the error term o(l), (|3.7p remains valid for the ranges of x\ 
and £2,3 up to o(e3) and o(es), respectively. 

In the remainder of this section, we prove Theorem 3.1. By expanding h(x, a) into finite Taylor 
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sum with the reminder term, we have 

x\ = -x\ + ^2 aij(a)xiXj + ^ a ijk (a)xiXjX k + 0(4), 
x 2 = a(a)x 2 - b(a)x 3 + ^ by (a) X{Xj -\- 

bi jk (a)xiXjX k + 0(4), (3.10) 
x 3 = b(a)x 2 + a(a)x 3 + V] Cij(a)xiXj + V] ^(a) + 0(4), 



where 



3 3 

i=l j=i 
3 3 3 

^j£ryjfc(a)xiXjXfc = > y aij k (a)xiXjX k , a£{a,b,c}. 

i=l j=i fc=j 

To study system of equations (|3.10p in a small neighborhood of the origin, we rewrite it in cylindrical 
coordinates (|3.2|) and rescale variables 

xx = £e 2 and p = re. (3-H) 

In new coordinates, we have 

£ = -^ + r 2 T 1 (9) + 0(e), 

r = er 2 R 1 {e) + e 2 r(p + ^R 2 (e) + r 2 R 3 {e))+0(e 3 ), (3.12) 
= l3 + erL 1 {9) + 0{e 2 ) , 

where Ti,L\ and #123 are homogeneous trigonometric polynomials. 
Using 4> = P _1 as a new variable, from (|3.12p . we obtain 

-^ + r 2 Px(^ + eE^,r, </»,e), (3.13) 



^ = eQ^r 2 + e 2 r ( M + £Q 2 (0) + r 2 Q 3 (</>)) + e 3 * (£, r, e) , (3.14) 

where Qi,2,3(<^) and Px((j>) are trigonometric polynomials of period u = 2irf3~ l . Functions E and 
Vl/ are bounded and continuous in 1/ x ]R x [0, 1]. Here D' stands for the domain of the rescaled 
variable (£,r) when the original variable x £ D C IR 3 : 

£>' = {(£,r) : \C\ < 4M 2 , < r < 2M} . (3.15) 

Similarly, we define 

D' Q = {(£,r): \£\<M 2 , < r < M} . (3.16) 

To determine the slow manifold for (|3,10p . we introduce the following notation. By pi and Px(<fi) 
we denote the mean value and the oscillating part of the trigonometric polynomial Pi(4>) in fj3. 13j) . 
respectively: 

px = A f * Pi(0)# and Px{4>) = Px{4>)-px, (3.17) 
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and define 



6>(0)=Pi+Pi(0) and Pl ^) = ^M—^l. (3.18) 

1 + Ap z 



We first prove an auxiliary lemma, which characterizes the trajectories of ()3.13j) and (|3.14p in 
D'. In the following two lemmas, we will keep track of how the remainder terms in the asymptotic 
expansions depend on the size of D'. This will be used later to determine the domain of validity 
for the asymptotic approximation of the slow manifold. Below k = O (Me) means that there exist 
positive constants eo and C (independent of e) such that \n\ < CMe for e € [0, eo]. In the remainder 
of this section, all estimates are valid for sufficiently large M > independent of e (as stated in the 
Theorem 3.1), but also allow the possibility that M > 1 grows as e — > 0. 

Lemma 3.1 Suppose (£,r) remains in D' for 4> e [0o,0]. Then 

m = £o(0)r 2 (0) + e -(*-*0 (£(0 O ) - £ o (0o)r 2 (0o)} + O (M 3 e) , € [fa, 4>] , (3.19) 

WHERE £o(<fi) IS DEFINED IN (13.181) . 

Proof: Using ([3.17)) . we integrate (|3.13j) over [4>o,4>] 

e*£(0) = e*°£(0o)+Pi f e s r 2 (s)ds+ f ' e s r 2 (s)P 1 (s)ds + e I* e s E (£(s), r(s), s, e) ds. (3.20) 
Using integration by parts, we have 

e s r 2 (s)ds = eV(0) - e^°r 2 (0 o ) - 2 f e s r(s)r(s)ds. 

J (p 

From (l3~T4l) and (f3T5j) . we find that for (£,r) G D' 

e s r(s)r(s)ds < Cie^°M 3 e, 
where Ci > is independent of e. Therefore, 

> 

eV(s)ds = eV(0) - e*V 2 (0 o ) + e^ O (M 3 e) . (3.21) 

Similarly, using integration by parts in the second integral on the right hand side of (|3.20p . we 
obtain 

f\ s r 2 {s)P l {s)ds = eV(<A)A(0)-eK 2 (0 o )A(0o) 

-2 f^e s r(s)f(s)Pi(s)ds- f ' e s r 2 (s)P[(s)ds. (3.22) 
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Applying integration by parts in the last integral on the right hand side of (|3.22p . we have 
e s r 2 (s)P{(s)ds = eV(^^)-e*r 2 (^)P[(^) 



e s r{s)r{s)P[(s)ds - I* e s r 2 {s)P[' (s)ds. (3.23) 

By direct verification, one finds that P" = (l + 4/3 2 ) P\. The integrals involving r on the right hand 
sides of (1332]) and (13331) are bounded by e^ O(M 3 e) since, by (|3T4j) . r = 0{eM 2 ). Using these 
observations, from (|3.22p and (|3.23p we obtain 

e s r 2 (s)P 1 (s)ds = eV(0)pi(0) - e*Y 2 (0o);pi(0o) + e^ O(M 3 e), (3.24) 

<t>o 

where pi(0) is given by (pU8]) . The combination of (|330jh (^2T) . and (13341) yields (j3T9]) . 
□ 

The following lemma shows that the trajectories of (|3.13|) and (|3.14p . which stay in D' for 
sufficiently long time enter a small neighborhood of S no later than in time 0(|lne|) and remain in 
this neighborhood as long as they stay in D' (see A3. 151) ) . 

Lemma 3.2 Let (£(0), r (0)) denote a trajectory of (13.1 3D and (13.141) . which stays in 
D' for cf> G [4> , 4>], 4> > 0i, 0i = </>o + | In e | . Then 

a0)=^o(0)r 2 (0) + O(M 3 e), (3.25) 

FOR > 01 AND AS LONG AS (£(0),r(0)) € £>'. 

Proof: Denote 0i = 0o + |lne|. By plugging in = 0i into (|3.19p . we have 

|^(0i)-eo(0i)r 2 (0i)| <C7 2 e, (3.26) 
for some C2 > independent of e > 0. Using Lemma 3.1 again with 0o := 0i, we obtain 

£(</>) = ^o(0)r 2 (0) + e~^^ (e(0i) - ^o(0i)r 2 (0i)} + 0(M 3 e), > X . 

The expression in the curly brackets is 0(e) by (|3.26p . 
□ 

Remark 3.2 Lemmas 3.1 and 3.2 show that the trajectories of (|3.13|) and (|3.14|) converge to an 
exponentially stable manifold S, whose leading order approximation is given in the definition of Sq. 
The method, which we used in Lemmas 3.1 and 3.2 to obtain the leading order approximation of 
the slow manifold, can be extended to calculate the higher order terms in the expansion for S. 

Having shown that the trajectories approach an O(e) neighborhood of So in time O (|lne|), next 
we reduce the dynamics of (|3.13p and (|3.14p to the slow manifold. For this, we define 

A0)=(^ + «/i(0)) 2 , (3.27) 
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where q'i{4>) = Qi{<fi) and Qi((f>) is a trigonometric polynomial on the right hand side of (|3.14p . For 
the sake of definiteness, we choose q'i(4>) such that 



/ qi(s)ds = 0. 
Jo 

The following lemma provides the desired reduction. 



Lemma 3.3 For <p > and as long as (£,r) G D', (I, 0) SATISFIES THE FOLLOWING SYSTEM 
OF EQUATIONS 

/ = -2e 2 (^/ + 7 + Q((/ ) )) + 0(M 3 e 3 ), (3.28) 
<j> = l + 0(Me), (3.29) 

WHERE Q(0) IS A TRIGONOMETRIC POLYNOMIAL WITH PERIOD CJ = 27T/? -1 AND ZERO MEAN 

f£ Q(4>)d<t> = 0. The expression for 7 is given in ( 13. 6j) . 
Proof: The change of variables 

r( ^) = t7j7\ 77\ ' where = ^i^)' ( 3 - 30 ) 

J(0) - egi (0) 

in Equation f)3. 14|) yields 

j = ^j- ((// + f Q 2 (0)) J 2 + Q 3 (0)) + 0(M 2 e 3 ). (3.31) 
After another change of variables, I = J 2 , we obtain 

g = -26 2 ((// + £Q 2 (<£)) / + Q-M) + 0(M 3 6 3 ), (3.32) 
By Lemma 3.2, for 4> > (pi 

m=U4>)r 2 {(t>)+0{Mh). (3.33) 
By plugging in (|3.33p into (|3.32[) . we obtain 

dl _ o , _ t , „ , 33, 



-2e z (»I + 6(0)^2(0) + Qs(<l>)) + 0(MV). (3.34) 
Equation (|3.28p follows by rewriting (|3.34p 

g = -2e 2 ( M I + 7 + Q(0)) + 0(M 3 e 3 ), (3.35) 

where ^ ^ 

7 = - / (6(0)Q2(0) + Qs(0)) # and / Q(<£)# = 0, (3.36) 
<*> Jo Jo 

and £o(^0 is given in (|3.18p . Equation (|3.29p follows from the last equation in (|3.12p and the 

definition of cf). 

□ 
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The statements in Theorem 3.1 can now be deduced from Equations (|3.28[) and (|3.29[) . We only 
need to show that a trajectory with initial condition in D' remains in D' for times longer than 
0(|lne|). This follows from the fact that 1 = (e 2 ) in D' . Consequently, it takes time 0(e~ 2 ) 
for / to undergo O(l) change necessary for leaving D' from D' . We omit any further details. In 
conclusion, we note that the domain of validity of the asymptotic analysis of this section extends 
much further beyond D' . Indeed, from f)3.25j) we observe that the remainder term tends to with 

e — > provided M = o ( e~ ). This means that the expansions for r and £ can be controlled in 



regions of size up to O (e 3 J and O (e a j respectively. For the original variables p, x\, these 
estimates translate into O ( e5 ) an d O respectively. 

We end this section by deriving a useful estimate for the time of flight of trajectories passing 
close to the saddle- focus. For this, consider an initial value problem for system of equations (|3. 13|) 
and (|3.14p . Suppose that the initial condition implies that I ((fro) = Iq > 0. We would like to know 
how long it takes for I to reach a given value 1 = 1 (0o + A), A > 0. We assume that Iq and I are 
sufficiently separated, e.g., 21 < Iq. To estimate A, note that 

A = Ai + A 2 , (3.37) 

where Ai = O ( | In e | ) is time necessary to reach an 0(e) neighborhood of So from Jo- Denote 

h :=/(0i) = / + O(e 2 |lne|) , where fa = O + Ai. (3.38) 

The second term on the right hand side of (|3.37p can be estimated by integrating (|3.28p over 
[</>!, 01 + A 2 ] : 

7 = I ie - 2aA * + 1(1- e~ 2aA2 ) - 2e 2 e - 2a ( A2 " s )Q (0 O + s) ds + O (e 3 ) . (3.39) 
M Jo 

The integral on the right hand side of (|3,39p is bounded uniformly for A 2 > 0. This observation 
combined with (|3.38[) implies 

I = I oe - 2aA > -1(1- e~ 2aA2 ) + O (e 2 ) . (3.40) 
A* 

Note that the contribution of 0(e 2 |lne|) term in (|3.38p to (|3.40p is negligible provided A 2 is 
sufficiently large. From (|3.40p . we obtain the desired estimate 



4 The oscillations 

In the present section, we use the local analysis of Section 3 to study certain oscillatory patterns 
arising in the model of solid fuel combustion (|2.ip - (|2.3p . By plugging in the values of the parameters 
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of (|2.ip - (|2,3p into the expression for 7 (|3.6p . we find that j(p) is a quadratic function with two 
zeros at p\ « 0.34 and P2 ~ 2.73. These values are in a good agreement with those obtained by 
numerical bifurcation analysis in [T3]. The quadratic character of j(p) is explained by the fact that 
7 is determined by the second order terms in the Taylor expansions of ^2,3(0, 0). We concentrate 
on the parameter region around p = P2, where 7(7?) changes its sign. For small positive a, there are 
two dynamical regimes: for values of p lying to the left and to the right of some 0{e) neighborhood 
of P2 corresponding to the supercritical and subcritical AHBs. The transition region between these 
two parameter regimes adds to the repertoire of qualitatively distinct dynamical behaviors. Below, 
we study these three cases in more detail. 

4.1 The supercritcal AHB: geometry of the periodic orbit 

It follows from Theorem 3.1 that for "f(p) < and sufficiently small a > 0, (|2.5p has a stable limit 
cycle O a , whose leading order approximation is given by 



(see Figure [6^). Moreover, the analysis of Section 3 shows that all trajectories starting from a suffi- 
ciently small neighborhood of the origin and not belonging to W s (0), converge to O a . The leading 
order approximation of O a in (|4.ip reveals a remarkable property of the oscillations generated by 
the limit cycle born from the supercritical AHB: the frequency of oscillations in x\ is twice as large 
as that of the oscillations in £2,3. To explain the frequency doubling effect, we recall that 



Therefore, (|4.ip implies that, unless A = 0, the frequency of oscillations in x\ is /Svr" 1 , while that 
of oscillations in £2,3 is J^. The latter coincides with the frequency of the bifurcating periodic 
solution. Below, we complement the analytical explanation of the frequency doubling with the 
geometric interpretation. 

The oscillations in 2:2,3 can be well understood using the topological normal form of the AHB 
[1, 20, 28j. Indeed, since at the bifurcation the center manifold at the origin is tangent to the X2 — £3 
plane, a standard treatment of the AHB using the center manifold reduction [7J [20J shows that the 
projection of the bifurcating limit cycle onto %2 — £3 plane to leading order is a circle (Figure 
[6b) and the projection of the vector field is given by the equation of the angular variable (|4.2p . 
Therefore, the oscillations in #2,3 are approximately harmonic with the period equal approximately 
to 2irf3~ 1 . The topological normal from, however, does not describe the oscillations in x%. For this, 
one needs to take into account the geometry of the bifurcating periodic orbit as a curve in IR 3 . The 
geometry of O a is fully determined by the two geometric invariants: the curvature, k(#), and the 
torsion, k(9). For the purposes of the present discussion, we need only the latter, but we compute 
both invariants for completeness. After some algebra, the parametric equation for the periodic 




(4.1) 



9 = /3 + 0(e). 



(4.2) 
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Figure 6: A periodic orbit born from the supercritical AHB (a), and its projections onto X2 — x-^(h) 
and x\ — X2 planes (c). Note that the periodic orbit shown in (a) is not planar. This accounts for 
the presence of the self-intersection in the its projection in (c) (see text for details). 



orbit P~T|) yields 



k(0) 
k(0) 



= + (3 cos 49 + 5) + h.o.t 

\x\ 
(A, x, x) 



"7 



QA sin 26 



x, x 



a l + 2,4 2 (3cos46> + 5) 



+ h.o.t, 6»€[0,2tt). 



(4.3) 



The geometry of the leading order approximation of the periodic orbit is determined by the three 
parameters a, 7, and A. The former two parameters are the same as the parameters in the 
topological normal form of the nondegenerate AHB [28J; while the latter captures the geometry 
of the slow manifold (or unstable manifold) near the origin. From the geometrical viewpoint the 
bifurcation is degenerate if either 7 or A is equal to zero. The latter condition holds if and only if 
the slow manifold is either a plane or a circular paraboloid near the origin. In this case, Equation 
(|4.3p implies that (to leading order) the bifurcating orbit lies in a plane. Generically, O a is not 
planar (see Figured^). The fact that the periodic orbit is generically not planar combined with 
the symmetry of the orbit about the xi— axis implies that its projection onto any plane containing 
27— axis has to have a self-intersection (see Figure [6b). From Figure [6)3, it is clear that as the phase 
point goes around O a once, x\ has to trace its range at least twice. Therefore, the frequency of 
oscillations in 27 has to be at least twice as high as that of the oscillations in 2:2,3. The above 
discussion implies that the frequency doubling of oscillations in x\ is a generic geometric property 
of the AHB. 



4.2 The subcritcal AHB: multimodal oscillations 

If the AHB is sub critical (7 > 0) the loss of stability of the equilibrium at the origin of (|2.5[) results 
in the creation of multimodal trajectories which spend a considerable amount of time near a weakly 
unstable equilibrium. To describe the resultant dynamics we give the following definition. 
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Definition 4.1 We say that a trajectory of (|2.5|) undergoes multimodal oscillations 

FOR t > IF THERE EXIST POSITIVE CONSTANTS T\ AND r 2 INDEPENDENT OF £ AND AN UN- 
BOUNDED SEQUENCE OF TIMES 

< t\ < t 2 < ■ ■ ■ , lim ti = oo, 

i— >oo 

SUCH THAT 

a. \x{t 2i )\ < ne 2 & \x{t 2 i-i)\ > r 2 , i = 1, 2, . . . , 

b. (3t' > h : \x(t')\ < rie 2 ) 3i £ IN : \x(t)\ < ne 2 ,t £ [min(t', t 2i ), max(i',i 2 i)] Vi' G 

(*2i-l) *2i+l) • 

The time intervals n = t 2 i+i — ^2i-i , * = 1,2,... are called ISIs. 

The proximity of (|2.5p to the AHB alone is clearly not sufficient to account for the appearance of 
the multimodal oscillations in (|2.5j) . Below, we formulate two additional assumptions on the vector 
field outside of the small neighborhood of the equilibrium, which are relevant to (|2.ip - (|2.3p . Under 
these conditions, we show that the subcritical AHB results in sustained multimodal oscillations. In 
addition, we determine the asymptotics of the ISIs for positive a = 0(e 2 ). For (|2.5p to generate 
multimodal oscillations, it is necessary that the trajectories leaving an 0(e 2 ) neighborhood of the 
origin reenter it after some interval of time. Therefore, the vector field in O(l) neighborhood of 
the origin must provide a return mechanism. Our second assumption on the global vector field 
of (|2.5p is that the trajectories approaching the origin along a ID stable manifold, W s (0), are 
subject to a strong contraction toward W s (0), i.e., in a small neighborhood of an 0(1) segment of 
W s (0), the projection of the vector field onto a plane transversal to W s {0) is sufficiently stronger 
than that along W s {0) (Figure [2h.). This guarantees that the trajectories entering such region 
of strong contraction approach W s {0) very closely and follow it to an 0(e 2 ) neighborhood of the 
unstable equilibrium, from where they are propelled away along the unstable manifold, W u (0). 
Due to the proximity of (|2.5p to the AHB, the motion away from the origin is very slow. This 
results in the pronounced intervals of the small amplitude oscillations. Below, we summarize these 
observations into two formal assumptions on the global vector field (Gl) and (G2). For this, we 
first need to introduce some auxiliary notation. For analytical convenience, we assume that in an 
O(l) neighborhood of the origin, W s (0) can be and has been straightened by a smooth change of 
coordinates. More specifically, the nonlinear terms in (|2.5p satisfy 

h-2,3 (xi,0, 0, a) = 0, w(xi,a) = —x\ + hi (xi, 0, 0, a) > 0, x\ E [di, 0), (4-4) 

for some d\ < 0. Note that (<ii,0,0) is assumed to be sufficiently far away from the origin (see 
Figure [8]). To describe the mechanism of return, we introduce two crossections: 

= {-c 1 5 2 ,c 1 5 2 ]x D s and E~ = {di} X D Sl , (4.5) 

where D$ C IR 2 denotes a disk of radius 5 centered at the origin: 

D s = {y E IR 2 : \y\ < 5} . (4.6) 
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Positive constants 5 and 5\ are sufficiently small (see Figure [2k) and c\ > is chosen so that S + 
intersects the slow manifold S transversally. In addition, we require that 5 = o(e 2 / 3 ) to guarantee 
that S + belongs to the region of validity of the local analysis of Section 3 (see Remark 3.1c). Let 
Xo G S + and consider a trajectory of (|2,5p starting from xo- We assume that for sufficiently small 
e > and a = 0(e 2 ), every such trajectory intersects S~ from the left. Denote the point of the 
first intersection by Q(xq) G £~. We assume that 

(Gl) the first return map Q : S~ — > S + depends smoothly on e and min E€ Q( S +) \x\ > C > 0. 
To measure the rate of contraction toward W s (0), we consider a 2 x 2 matrix 

/ 9/2 9/2 \ 

A(*i,a)= || • (4-7) 

V dx 2 dx 3 J ( Xli0 ,0,a) 

Let Ai(xi,a) < \2{xi,a) denote the eigenvalues of the symmetric matrix 

A s ( Xl , a) = - (A (xi, a) + A T (x u a)) . 

Denote X(xi,a) = —\i(xi,a) and \(xi,a) = —\2{xi,a). We assume that for sufficiently small 
e > and a = 0(e 2 ), 

(G2) 

^-A(xi,a)<0, x 1 €[d 1 ,0]] (4.8) 

OX\ 

3d 2 € (di,0) : X(xi,a) > 0, x 1 e[d 1 ,d 2 ] & min X ^ a \ =0(|lne|), (4.9) 

xie[di,d2]w(xi,a) 

(G3) 

max = o(|i ne n 

xiG[rfi,o) w(xi, a) 

Under these conditions, we have 
Theorem 4.1 Let 7 > 0, conditions in (14.41) . (Gl) and (G2) hold. Then for sufficiently 

SMALL e > AND a = O (e 2 ) , A TRAJECTORY OF (12. 5|) WITH INITIAL CONDITION FROM AN 
0(e 2 ) NEIGHBORHOOD OF THE ORIGIN AND NOT BELONGING TO W s (0) UNDERGOES MULTIMODAL 
OSCILLATIONS. THE ISIS ARE UNIFORMLY BOUNDED FROM BELOW 

n>T' = ^\n(l + ^jC-) , i = 2,3, (4.10) 



2a ~ V ' 7e 4 

IF, IN ADDITION, (G3) HOLDS THEN THE ISIS SATISFY TWO-SIDED BOUNDS 

t- < n < t+ = -L In ( 1 + -^—C + V 4 = 2,3,..., (4.11) 



2a V 7 g4+x 

for some x > 0. Positive constants C 1 * 1 do not depend on e, a, and 7. 
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Figure 7: The interpretation of the bifurcation scenario arising in (|2.ip - (|2.3p for fixed positive 
value of a = 0(e 2 ) and varying 7 near using a one-parameter family of ID maps. For negative 
values of 7 = 0(1), the map in (|4.17j) has a stable fixed point (upper graph), which corresponds 
to a stable limit cycle born from the supercritical AHB. For increasing values of 7 the fixed point 
looses stability via a period-doubling bifurcation (lower graph), which initiates a period-doubling 
cascade leading to the formation of the chaotic attractor. 



Remark 4.1 



(a) The principal assumptions on the global vector field are formulated in (Gl) and (|4.9p . The 

condition in (|4.8p makes the derivation of certain estimates in the proof of the theorem easier 
and is used for analytical convenience. Likewise, (G3) is not essential for the proposed 
mechanism. However, without this condition obtaining the two-sided estimates for the ISIs 
requires an additional argument in the proof. Condition (G3) means that the two eigenvalues 
of A have the same order of magnitude. This condition is not restrictive. 

(b) For fixed 7 > 0, the estimates in (j4.11j) and (|4.10p can be rewrittten as 

Mi + g-q) < n < m(i + c + q) (4 12) 

2a ~ ~ 2a K J 

with constants independent from a. Similarly, for fixed a and varying 7 > inequalities 
in P~TT|) and (jlTTUj) imply 

0- - ^ln 7 < n < 0+ - -Lln7, (4.13) 
where constants = O ( |ln e| ) do not depend on 7 > 0. 



4.3 The transition from subcritical to supercritical AHB 

The numerical experiments presented in Section 2 show that the transition from subcritical to su- 
percritical AHB contains a distinct bifurcation scenario involving the formation of chaotic attractor 
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via a period-doubling cascade. The analytical explanation of these phenomena is outside the scope 
of the present paper. Below we comment on the difficulties arising in the analytical treatment of 
this problem. In the present section, we use a combination of the analytic and numerical techniques 
to elucidate the origins of the complex dynamics arising in the parameter regime near the border 
between regions of sub- and supercritical AHB. The principal features of this bifurcation scenario 
are summarized in Figures 0] and [5j In these numerical experiments, we kept a fixed at a small 
positive value and varied 7. By taking progressively smaller values of 7 > 0, one first observes 
that the multimodal patterns exhibit an increase in the ISIs (Figure H^,b; see also Figure [HJa). We 
consider these oscillatory patterns regular (even if they are not periodic), because the timings of 
the spikes remain within narrow bounds in accord with (|4.1ip . For smaller values of 7 > 0, the 
oscillatory patterns become irregular and are characterized by long intervals of oscillations of small 
and intermediate amplitudes between successive spikes (Figure H]c,d). As 7 becomes negative, the 
trajectories lose spikes and consist of irregular oscillations (Figure [5] a, b). Further decrease of 7 
leads the system through the reverse cascade of period-doubling bifurcations (Figure [5] c-h). The 
period-doubling cascade terminates with the creation of the limit cycle, which can be followed to a 
nondegenerate supercritical AHB by letting a — > (Figure [5] g,h) . To account for the bifurcation 
scenario described above, we construct the first return map. Since near the origin the trajectories 
spend most of the time in the vicinity of the 2D slow manifold, the first return map is effectively 
one-dimensional. The distinct unimodal structure of the ID first return map affords a lucid geomet- 
ric interpretation for the bifurcation scenario in the 3D systems of differential equations near the 
transition from sub- to supercritical AHB. Specifically, we show that the mechanism for generating 
complex dynamics in the continuous system in this parameter regime is the same as in the classical 
scenario of the period-doubling transition to chaos in the one-parameter families of unimodal maps 
[20]. 

We next turn to the derivation of the first return map. We start with extending the asymptotic 
analysis of Section 3 to cover the case of 7 = o(l). This requires computing additional terms on the 
right hand side of the reduced equation (|3.28j) . because already for [7] = 0(e), the term involving 
7 on the right hand side of (13.281) is comparable with the 0(e 3 ) remainder term. We then use the 
more accurate approximation for the equation for / to compute a ID mapping: 

G : 1(0) ^ I (0 + 27T/T 1 ) , (4.14) 

which describes how / changes after one cycle of oscillations. As noted above, the reduced equation 
(I3.28j) is not suitable for analyzing the case of small [7! and one needs to include more terms in 
the expansion on the right hand side of (|3.28p . In analogy with the topological normal form for 
the degenerate AHB [28l [19], we expect that terms up to 0(e 5 )) are needed in (13.28H to resolve 
the dynamics for small |7|. This can be achieved by a straightforward albeit tedious calculation. 
Below we describe the formal procedure of obtaining the required expansions. The justification of 
these expansions can be done in complete analogy with the analysis of Section 3. First, we compute 
terms on the right hand sides of (|3.13p and (|3,14p up to 0(e 3 ) and 0(e 6 ) respectively. Then we 
look for a solution of A3. 13H in the following form: 

C(<t>) = £o(r, 0) + e£i(r, 0) + e 2 £ 2 (r, <j>) + 0(e 3 ). (4.15) 

By taking an initial condition from Dq and using Anzats (I4.15P in A3. 131) . we recover (,o(r,(j)) (see 
(|3.18p ) and find the next two terms in the expansion of £(r, (f>): £i(r,(ft) and £2(^5 </>)• Thus, we 
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Figure 8: The plots of first return map for X (|4.19p computed for (a) a = 0.0445, 7 = —0.0486, 
(b) a = 0.0445, 7 = -0.0126, (c) a = 0.0445, 7 = 0.0290, (d) a = 0.12, 7 = -1.0741. The plots 
in (a-c) illustrate the bifurcations in the family of the first return maps for the values of 7 near 
zero. The map shown in (d) corresponds to the irregular oscillations following the period-doubling 
cascade for increasing values of a in the case of supercritical AHB (see Remark 4.2). 
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obtain the approximation of the slow manifold with the accuracy 0(e 3 ). Next, we plug in (|4.15p 
into (|3.14p and collect terms multiplying equal powers of e to obtain 

r = eR l (r, & + ■■■ + e 5 R 5 (r, <j>) + 0(e 6 ), (4.16) 

where Rk(r,4>) are 27T/3 -1 periodic functions of the second argument. Finally, by rewriting (|4.16p 
in terms of I (see (3.25)), integrating it over one period of oscillations, u, and disregarding 0(e 5 ) 
terms, we obtain a map for the change in I after one cycle of oscillations 

I n+1 = G(I n ), G(I) = I — 2e 2 uo (^ l I + 1+ e -^ iUJ = 27T/T 1 . (4.17) 

Except for the last term in the definition of G, the map in (|4,17p follows from the reduced system 
(|3.28p and (|3.29p . The last term is only needed if the value of | 1 does not exceed 0(e). Our 
calculations show that the second Lyapunov coefficient, c, is negative for the values of parameters 
used in (|2.ip - (l2.3p . Therefore, for sufficiently small e > 0, (|4.17p defines a unimodal map. Away 
from an 0(e 2 ) neighborhood of 0, the graph of / = G(I) is almost linear with a weakly attracting 
slope 1 — 2iope 2 for p > (Figure [7|). For p > 0, G has a unique fixed point: 

^+je 2 + 0(6 3 ), 7 <0, 

fe + 0(e 2 ), 7 = 0, (4.18) 
: e 2 + 0(e 3 ), 7 >0. 
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For negative values of 7 = 0(1), I = lis a stable fixed point of G, which corresponds to a stable limit 
cycle born from a supercritical AHB (see Figure [5] g,h). It follows from (|4.17p that for increasing 
values of 7, the graph of the map moves down. Consequently, the fixed point moves to the left 
and eventually looses stability via a period-doubling bifurcation (Figure [5] e,f). Further increase 
in 7, yields more period-doubling bifurcations (Figure [5] c,d). Given the unimodal character of G, 
one expects that this sequence of period-doubling bifurcations eventually leads to the formation 
of a chaotic attractor (Figure [5] a,b). The unimodal character of the map combined with the 
manner of its dependence on 7 suggest a clear geometric mechanism for the formation of the 
chaotic attractor and the subsequent period-doubling cascade arising near the transition from sub- 
to supercritical AHB (Figure [5h,-h) . We now outline the limitations of the asymptotic analysis of 
this section and difficulties arising in the justification of (|4.17p . According to (|4.18p already at the 
moment of the first period-doubling bifurcation, fixed point / belongs to an 0(e 2 ) neighborhood 
of 0. In this neighborhood, p = 0(1) (see (|3.27p ) and, therefore, I lies outside of the region of 
validity of the asymptotic analysis. In this region, (|4.17p may only be considered as a formal 
asymptotic expression. A rigorous justification of the from of the map in the boundary layer 
meets substantial analytical difficulties: it requires introducing an additional set of intermediate 
asymptotic expansions and matching them with those obtained in Section 3. We do not address 
this problem in the present work. Below, we resort to using numerical techniques to verify the 
principal features of the first return map suggested by the asymptotic analysis: the unimodality 
of G and its dependence on 7. The numerics confirmed our predictions about the form of the first 
return map and it also revealed certain additional features. 

For numerical construction of the first-return map, we fix angle 9 E [0, 2ir) and define a 
crossection in cylindrical coordinates: £ = {(p, £,6>) : 6 = 8}. Let (p(t),^(t),9(t)) be a trajec- 
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tory of (|2.5p and < t\ < ti < t% < • ■ ■ < tk < • • • denote a sequence of times, at which 
(p(ifc), 0(i&)) G X, k = 1, 2, 3, .... Then, we define 



£ : Z(ifc) i-» and J(i) = 




(4.19) 



By a suitable choice of 0, we can achieve that a trajectory of (12. 5p . intersects S once during one 
cycle of oscillations. Therefore, for a trajectory lying in a 2D slow manifold, (|4.19p defines a 1.D 
first-return map. Note that maps (|4.17p and (14.19P are related via rescaling of the coefficients, since 
to leading order X = e~ 2 I. Therefore, the graphs of G(I) and Q(I) are similar in the domain, where 
(|4.17p is valid. We computed the first return maps for fixed a = 0.0445 and for several values of 
7. The representative plots are shown in Figure [8] a-c. The numerically computed maps confirmed 
our expectations about the graph of Q in the boundary layer near 0: in this region, the map is 
decreasing with a strongly expanding slope. Away from 0, the graph of Q contains an almost linear 
branch, whose slope is very close to 1 — 2atu, as predicted by (|4.17p (see lower branches in Figure 
[Hb, c). These numerics also confirm that under the variation of 7, the graph of the first return map 
is translated in vertical direction (Figure [8] a-c). Therefore, the family of the first return maps 
possesses two principal ingredients (the unimodality and the additive dependence on 7), which are 
necessary for the qualitative explanation of the bifurcation scenario given in Figure [7J In addition, 
our numerical experiments reveal a new feature of the first return map: for small I7I, the graph 
of the map has another almost linear branch in the outer region away from the origin (see Figure 
(lb)- This branch of the graph of the map also has an attracting slope. The presence of this branch 
in the first return map indicates the existence of another (branch of the) slow manifold different 
from that described in Section 3. A possible explanation for the appearance of the upper branch 
in the first return map is due to the unstable manifold of the periodic orbit. For small I7I, the 
first-return map (|4.19p is multivalued in the outer region. However, since both the upper and 
the lower branches have positive slopes less than 1, the qualitative dynamics for the map shown 
in Figure [8b does not depend on the exact mechanism for the selection between the branches in 
(|4,19p . Although explaining the multivaluedness of the first return maps shown in Figures [8b, c 
presents an interesting problem, it is not critical for the qualitative explanation of the bifurcation 
scenario arising during the transition from the subcritical to supercritical AHB. The latter was the 
main goal of the present subsection. 

Remark 4.2 In addition, to the region in the parameter space containing the border between the 
regions of sub- and supercritical AHB, there is another parameter regime resulting in the complex 
dynamics. In [13], it was shown numerically that the limit cycle born from the supercritical AHB 
in (|2.5p undergoes a period-doubling cascade leading to the formation of the chaotic attractor for 
increasing values of a. The first period-doubling bifurcation in this cascade is shown in Figure 
[Tb,d. This bifurcation scenario is consistent with the form of the first-return map constructed in 
this subsection. Indeed, it is easy to see from (I4.18P that for increasing values of fi the fixed point, /, 
moves to the left. Therefore, the explanation given above for the period-doubling cascade resulting 
from the variation of 7 near also applies to the case of increasing a and fixed 7 < (see Figure 

EH). 
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Figure 9: A schematic representation of the trajectory of (|2.ip - (|2.3p approaching a weakly unstable 
saddle-focus. The trajectory crossing E~ is subject to strong vector field transverse to W s (0) in 
Tig . This guarantees that it remains close to W s (0) until it enters an 0(e 2 ) neighborhood of the 
origin. From that moment and until the trajectory reaches S + its behavior is described by the 
local analysis of Section 3. 



5 The proof of Theorem 4.1 



In the present section, we show that under the assumptions of Theorem 4.1 the trajectories of (??) 
with initial conditions in S~ enter Dq. The local analysis in Section 3 describes the behavior of 
trajectories from the moment they reach Dq and until they leave D. In particular, it shows that 
the dynamics near the origin has two phases: the fast approach to the slow manifold and the slow 
drift away from the origin along the slow manifold. Upon leaving D, the trajectories are reinjected 
back to £~ by the return mechanism postulated in (Gl). This scenario implies that the system 
undergoes multimodal oscillations as stated in Theorem 4.1. 

We start with presenting several auxiliary estimates, which will be needed for the proof. By 
(|4.4j) . one can choose C\ > such that 

w(xi,a) > C\ \x\\ , x\ £ [<ii , 0] (5-1) 

and for sufficiently small a > 0. Next, we note that A(0, 0) = 0. This equation and (|4.8j) . by the 
Implicit Function Theorem, imply that there exists > x*(a) = 0(e 2 ) such that 



X(xi, a) 



<0, x\ € [di,x*(a)), 
> 0, x\ € {xl(a),0]. 
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Using (|4.8p and (|5.ip . we have 

A(xi,a) A(xJ,a) 
w(xi,a) w(xi,a) 

where w(xi,a) is defined in (|4,4I) , 

Let // > be such that (Gl) and (G2) hold for a € [0, /ie 2 ]. In the remainder of this section, 
unless stated otherwise, it is assumed that a G [0, fie 2 ] . To simplify notation, we will often omit 
the dependence of various functions on a. Using (|4.4p . we rewrite (|2.5p in the following form: 

ii = io (x x ) +yi<AiOi>y) + 2/2<fo (#1, y), (5-2) 
y = + ip(xi,y), (5.3) 

where to(xi) and A(x\) are defined in (14. 4p and (14. 7p respectively; y = (yi,y2) := (^25^3), and 

w(0) = 0, 1>(x x ,v) = O (|y| 2 ) , 0i,2(xi, y) = 0i,2(si) + 0(\y\) , x £ [di, 0] . (5.4) 

Let 

d 3 = -M 2 e 2 (5.5) 

denote the xi— coordinate on the left lateral boundary of Dq (see ()3.3[) ) and LT^ = [^1,^3] x .D^ 
(see Figured]). From (|2.5f) and (|4.4p . one can see that for 5\ > sufficiently small (independent of 
e), the right hand side of (|5.2p 



+ yi0i(^i,y) +y2<h(xi,y) > o, (xi,y) e n^. 

Therefore, in II ( 5 1 , (15. 2p and (15. 3p may be rewritten as follows 



dy(x 



dx 

where 



^=i(xi)l/ + fe,y), (5-6) 



i(x x ) = and ^( Xl ,y) = ^^(l+O (|y|))- 

w(^i,y) 



By (fSTTJ) and JO) , we have 



C I I 2 

Hxi,y) < f , , (xi,y)en 5l! (5.7) 

for some C2 > independent of e. 

To follow the trajectories from S~ to Dq, we introduce two regions: 

11^ = [di,d 2 ] x D 5l and II 2 2 = [c^cfe] x D s 2 (see Figured]). 

Recall that D$ l 2 denotes the disk of radius 5^2 (see ()4.6p ). Positive constant 5\ is the same as 
in (|4.5p and < 62 <C £1 will be specified later (Figured]). Recall that d% = —Me 2 denotes the 
value of the x\— coordinate on the the left lateral boundary of Dq (see (|3.3p ). By taking M > 
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large enough, we can arrange d 3 < x*(a) < for a E [0, /ie 2 ] . By (Gl), the vector field in L 7 ^ is 
sufficiently strong so that the trajectories entering ITj through S~ get into a narrow domain Il| 
and remain there until they reach Dq. The following lemma allows to control the trajectories in 
U 5l = [di,d 3 ] x Dg 1 . 



Lemma 5.1 Let [£i,£ 2 ] Q[di,dz], < 6 < Si, and 

max si 

xi€[ti,h] \y\ 



r — max sup yA(xi)y,yj < 0. (5i 



Then 

Wi)\ < 5 => \y( Xl )\ < 25e- v ^-^\ x x e [£ u £ 2 ], (5.9) 

PROVIDED 

4C 2 5 (1 - e -^-^)) 

Sn - L 5 - 10 

\£ 2 \ v 

Proof: We may assume that |y(xi)| / 0, x € [4, «Y], since, otherwise, by the uniqueness of solution 
of the initial value problem for (15.61) . y{x\) = 0, x\ € [^1,^2], and (15. 9h holds. From (|5.6p . (|5.7p . 
and (|5.8p . we have 



^^<-y|y(xi)| + ^(xi,|y|), ^( Xl ,\ y \) = ^f, x 1 e[£ 1 ,£ 2 ),\y\<\S 1 \. (5.11) 

where u is defined in (|5.8|) . and |y| denotes the Euclidean norm of y € IR 2 . Let y{x\) denote the 
solution of the initial value problem 

y' = -vy + i>{x 1 ,y),y{£ l ) = 5. (5.12) 

It is sufficient to show that (j5.9l) holds for y(xi). We represent y(xi) as the limit of the sequence 
of successive approximations 

y(™+!)( Xl ) = $ e -«(*i-<i) + jT* e- v{xi ~ s) ^ (a, y (n) (s)) da, n = 0,1,2,... (5.13) 

and y(°\x\) = 5e~ v ( Xl ~ il \ We use induction to show 

/ {n) (zi)| <25e-^ xi - £l \ xe[h,£ 2 ], n = 0,1,2,.... (5.14) 

Inequality (|5.14p holds for n = 0. We show that (|5.14p T ,— ^ =4> {5H3Dn=fe+i- Using the definition of 
'tp in (|5.1ip and (|5.14p n -fr. we have 

^e-^i-^a.j^sjds < ^ /" X1 e-^i— ) y( fc )( s ) 2 rf s <^ H e - v ^- s hS 2 e~ 2v ^-^ds 

l 1*2 1 J£i |«2 1 

= ^-.i^c^-;'"-'-'), (5 , 5) 

Using (I5.15j) and (|5.10p . from (|5.14j> ^_ t. we obtain (15.14p „-i% + i . By induction, (I5.14j) holds. By 
taking n — * oo in (|5.14p . we have 

|y(xi)| < 25e- v ^-^\ x x € [4,4]- (5.16) 
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The statement in (|5.9p then follows from (|5.16p and the standard theory for differential inequalities 

[22]. 

□ 

In the following lemma, we determine the size of Tig . 
Lemma 5.2 There exists positive 5 2 = 0(e 2 ), such that the trajectories of (|5.2|) and 

(15. 3P ENTERING Il| FROM Ug REMAIN IN ILg UNTIL THEY REACH Dq (SEE FIGURE [9j) 

Proof: Denote 

— V2 = max sup (A(xi)y, y) . 

xi€[d,24i] \y\=l ^ ' 

Recall 

X(x*) = and 0(e 2 ) = d 3 < x* = 0{e 2 ). 
Therefore, by (|4.8p . A(ds) = 0{e 2 ) is positive. In addition, 



V2 = — max sup [A(xi)y,y) = min 



•sup (y | =1 



xie[d 2 ,<fe] |j/|=l ^ ' / zie[d 2 ,d 3 ] M^l) 

min^g^^] A(xi) _ A(d 3 ) _ 2 



max 



ciefo,^] max zie[d 2 ,d 3 ] w ( x i) 



and f2 > 0. Next we apply Lemma 5.1 with 5 := 62, v := V2, and £1,2 := cfe,3- Note that for small 
V2 the inequality (15.10P can be rewritten as 

AC252 {d 3 -d 2 + Q(v 2 2)) < L 
1^3 

Since ^3 = 0(e 2 ) (see (|5.5p ) and u 2 = 0(e 2 ), one can choose 82 = 0(e 2 ) so that (|5.17p holds. 
□ 

Having found the size of Il 2 2 , we now determine the rate of contraction in IT] sufficient to 
funnel the trajectories entering ILg through X - to Hg 2 - 



Lemma 5.3 The trajectories of (15. 2D and (15. 3[) entering Ug through S remain in 

Tig UNTIL THE> 

Proof: Denote 



Tig UNTIL THEY REACH ILg . 

—v\ = max sup (A(x\)y,y 

xi£[di,d,2] [2/1=1 ^ 

From (G2), one finds that v\ = 0(|lne|) is positive. Lemma 5.3 now follows from Lemma 5.1 with 
5 := 5\, v := v\, and t\ t 2 := di,2- Indeed, for v\ = O (|lne|), we have 

vi > -rry-. (5.18) 

I "2 1 

Inequality (|5.18p is sufficient for (|5.10p to hold. By Lemma 5.1, we have 

\y(d 2 )\<26 1 e~ v ^ d2 - dl \ (5.19) 



30 



With vi = (|lne|), by (|5.19p . we can achieve |y(d 2 )| < S 2 = 0(e 2 ). 
□ 

Lemmas 5.2 and 5.3 imply that the trajectories entering through £~ stay in IlJ |J II 2 - until 
they reach Dq. Moreover, the inequality in (Gl) guarantees that such trajectories are bounded 
away from W s {0). Thus, we can use the analysis of Section 3 to describe the evolution of trajectories 
from the moment they reach Dq until they leave D. By Remark 3.1c, this description extends to 
any region where X2,3 = o (e 2//3 ) , i-e., the trajectories can be controlled until they hit S + . This is 
followed by the return to IB according to (Gl), and the next cycle of the multimodal oscillations 
begins. The analysis of this section applies to any trajectory starting from Dq and not belonging 
to W s (0). 

It remains to estimate the ISIs. For this, we compute the time needed for the trajectory starting 
in an 0(e 2 ) neighborhood to return back to this neighborhood after making one global excursion. 
Since the time of flight of the trajectory outside a small neighborhood of the origin depends regularly 
on the control parameters, the duration of the very long ISIs is determined by the time spent in 
the neighborhood of the origin. To estimate the latter, we note that at the moment a trajectory 
enters Dq, we have 

\y(tk)\ <C 3 e 2 . (5.20) 

This follows from Lemma 5.3, since [y^) < 62 and 62 = O (e 2 ) . To obtain the lower bound on 
[3/(^3) I , recall that by (Gl), we have 

|tf(di)|>C>0. (5.21) 

As follows from (G3), the maximal rate of contraction in 11^ does not exceed O ( |ln e| ) in absolute 
value. This combined with (j5.21j) implies that 

|y(d 3 )|>C 4 e 2+ f (5.22) 

for some x > and C4 independent of e. We omit the proof of (|5.22p . because it is completely 
analogous to that of Lemma 5.1. 

Let t = t~ denote the moment of time when a trajectory of (12. 5p enters Dq from II 2 2 . After 
switching back to the original parametrization of y by time, t, we rewrite (|5.20p and (|5.22|) : 

C 3 e 2+ ^ < \y{t~)\ < C 4 e 2 . (5.23) 

For t > t~ the trajectory approaches and remains close to the slow manifold as long as \y\ = o(e 2 / 3 ) 
(see Remark 3.1c). Let | < j < 1 and denote 

t+ = mm{t > t~ : \y(t)\ = e j ) . (5.24) 

From (|5.23|) and (|5.24p . we have the following bounds for Iq := I(t~) and I\ := I(t + ): 

C & e~ 2 -* <I < C 7 e- 2 and h = 0(e 2 ^), ie(|,l). (5.25) 
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For these ranges of values of Iq and I\ , from (|3.4ip we have 



r m = t + - t~ = In ( ( 1 + £l ) (l + ^o(l) ) ) . (5.26) 



2//e 2 \ \ 7 / V 7 
The combination of (|5.25|) and (15. 26ft yields two-sided bounds for Tj n : 

rr < r in < r+, (5.27) 

where 

.\a(l + JLc~) and r+ = -^ In ( 1 + — £—C+ 



m 2//e 2 V 7£ 2 / m 2^e 2 V 7^ 2+x 

where positive constants C% can be chosen independent of e. Inequalities (|5.27p provide bounds for 
the time that a multimodal trajectory spends near the unstable equilibrium. On the other hand, 
the time of flight outside a small neighborhood of the origin, r out , depends regularly on e, by (G2). 
Therefore, by adjusting constants C% in (|5.27p if necessary, one can obtain uniform bounds on the 
ISIs r = r in + r out for sufficiently small e > 0, positive 7 and n from bounded intervals, as stated 
in Theorem 4.1. 



6 Discussion 



In the present paper, we investigated a mechanism for generation of multimodal oscillations in a 
class of systems of differential equations close to an AHB. Our analysis covers both cases of sub- 
and supercritical AHB. For the supercritical case, we identified a novel geometric feature of the 
bifurcating limit cycle, the frequency doubling effect. It turns out that generically in the normal 
system of coordinates the oscillations in one of the variables are twice faster than in the remaining 
two variables. Therefore, the leading order approximation of the limit cycle bifurcating from the 
supercritical AHB requires two harmonics. The asymptotic analysis of the present paper explains 
the frequency-doubling. In addition, we provide a complementary geometric interpretation to this 
counterintuitive effect. In particular, we showed that it is a consequence of the geometry of the 
limit cycle. The latter is not captured by the topological normal form of the AHB. The analysis of 
the multimodal oscillations arising from the subcritical AHB requires additional assumptions on the 
global behavior of trajectories. We identified two principal properties of the global vector field: the 
mechanism of return and the strong contraction property. In the presence of this global structure, 
the subcritical AHB produces sustained multimodal oscillations combining the small amplitude 
oscillations near the unstable equilibrium with large amplitude spikes. The resultant motion is 
recurrent in a weak sense: it may not be periodic but nevertheless the timings of the spikes possess 
certain regularity. We have shown that the ISIs have well-defined asymptotics near the AHB 
and comply to the two-sided bounds, which depend on the principal bifurcation parameters. Our 
estimates show that near the AHB, the ISIs can be extremely long and can change greatly under 
relatively small variation of the bifurcation parameters. The ability of the system to exhibit such 
extreme variability in the ISI duration is important in many applications, in particular, in the 
context of neuronal dynamics. Previous studies investigated different possible mechanisms for 
generating multimodal patterns with very long ISIs [31 \TT\ \12\ [30], [35] . For the finite dimensional 
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approximation of the model of solid- fuel combustion (|2.ip - (|2.3p . the main motivating example 
for our work, the proximity to the homoclinic bifurcation was suggested in [13] as a possible 
mechanism for prolonged ISIs. Our conclusions confirm the importance of the proximity to the 
homoclinic bifurcation for explaining the oscillatory patterns in (|2.ip - (|2.3p . The proximity to the 
homoclinic bifurcation is implicitly reflected in our assumptions on the global vector field. However, 
we emphasize the critical role of the AHB: the duration of the ISIs can be effectively controlled 
by the parameters associated with the AHB without changing the distance of the system to the 
homoclinic bifurcation. In all our numerical experiments, the system remained bounded away from 
the homoclinic bifurcation, nevertheless it exhibited patterns with very long ISIs whose duration 
was amenable to control. Our analysis extends the estimate for the ISIs obtained in [21J to a 
wide class of problems. It also emphasizes the proximity of the system to the border between 
sub- and supercritical AHB, as another factor in creating oscillatory patterns with long ISIs. We 
show that as this border is approached from the subcritical side the ISI grow logarithmically. This 
observation is important for explaining the oscillatory patterns generated by (12.1|) - (j2.3l) since in this 
model the AHB can change its type under the variation of the second control parameter p. This 
situation is not specific to the model of solid fuel combustion. Recent studies suggest that there is 
a class of neuronal models close to the AHB whose type may change with the values of parameters 
133]. Therefore, it is important to understand the principal features of the transition from sub- 
to supercritical AHB. We found that when the border between the regions of sub- and supercritical 
bifurcation is approached from the subcritical side (while the distance from the AHB remains 
fixed), the oscillations become chaotic. The regime of irregular oscillations is then followed by the 
reverse period doubling cascade. To understand the nature of this bifurcation scenario we used a 
combination of analytic and numerical techniques. Using the insights gained from the asymptotic 
analysis, we constructed a ID first-return map. The map provides a clear geometric interpretation 
for the bifurcation scenario near the transition from sub- to supercritical AHB. Our study suggests 
that the formation of the chaotic attractor via a period-doubling cascade is a universal feature 
of this transition. For example, the bifurcation scenarios reported for the Hodgkin-Huxley model 
in [11] are very similar to those studied in the present paper and are likely to share the same 
mechanism. 

Mixed-mode oscillations similar to those studied in the present paper, have been studied in 
for a class of the slow- fast systems in IR 3 near the AHB. Although, the work toward developing a 
complete mathematical theory for such oscillations is still in progress, the general mechanism for 
their generation and the bifurcation structure of the problem have been greatly elucidated recently 
[261 l30l l3T| 138] . The present paper shows the relation between the mechanisms for the mixed- 
mode oscillations in the model in [13] and for those in the slow-fast systems. The latter possess 
a well-defined structure of the global vector field due to the presence of the disparate timescales 
in the governing equations [24] . The analyses of the mixed- mode oscillations in [26} [30| l3T| 155] 
use in an essential way the relaxation structure of the problem. The model in |13] is an example 
of the mixed-mode generating system, which does not possess an explicit relaxation structure. In 
fact, it is hard to expect such structure in a system obtained via projecting an infinite-dimensional 
system onto a finite-dimensional subspace. In formulating the assumptions on the global vector- 
field (Gl) and (G2), we were looking for the minimal requirements on the system near an AHB 
that guarantee the existence of the mixed-mode solutions. Due to the lack of the information about 
the global vector field of (|2.ip - (|2.3p . it appears impossible to verify these conditions analytically. 
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However, the numerical simulations clearly show that system of equations (|2.ip - (|2.3p possesses 
the qualitative structure required by (Gl) and (G2) (Figure [2b). On the other hand, we expect 
that conditions (Gl) and (G2) should be possible to verify analytically for a wide class of slow- 
fast systems. Therefore, we believe that our results will be useful for understanding mixed-mode 
oscillations in such systems. In particular, it would be interesting to apply this approach to the 
modified Hodgkin-Huxley system |TTJ. The numerical results reported in [llj strongly suggest that 
the mechanism proposed in the present paper is responsible for the generation of the very slow 
rhythms and chaotic dynamics in the Hodgkin-Huxley model. 
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Appendix A. 

In this appendix, we list explicit expressions of various constants and trigonometric polynomials, 
which appear in the definitions in of the slow manifold, S, and the first Lyapunov coefficient, 7. 
All expressions are given in terms of the coefficients of the power expansions on the right hand side 
of (|3.10p . By Oijk we denote Oijfe(O), a € {a,b,c}. The following constants are used to define the 
leading order approximation of the slow manifold in (|3.4p and (|3.5p : 



«22 + 033 , a 23 + (°22 - 033) j n , 2/?a 2 3 + (a 2 2 ~ 033) 

a = , A = \ — — A —. -— r- , and if = arctan — -r . 

2 V 4(l + 4/3 2 ) a 23 + 2/3 (a 22 - a 33 ) 

The following trigonometric polynomials enter the right hand sides of (|3.13|) and (13.140 : 

Pi (4>) = a 22 cos 2 (3(f) + a 2 3 cos (3(f> sin (3(f) + a 33 srn2 P&i 

Q\{4>) = b 2 2 cos 3 /30 + (6 23 + c 22 ) cos 2 /30sin/30 + (633 + c 23 ) cos p<f> sin 2 /3<fi + c 33 sin 3 , 
Q2{4>) = bu cos 2 @4> + (613 + C12) cos P4> sin f}(j> + C13 sin 2 f3<f>, 
Q3(4>) = b 2 22 cos 4 13(f) + (6 22 3 + c 222 ) cos 3 (3<fism.f3<fi + (b 2 33 + Q223) cos 2 /3</>sin 2 
+(^333 + c 23 3) cos /30sin 3 /3c/> + c 333 sin 4 13(f), 



Functions Q2,z{Q) = Q2,3(P l G) are used in the calculation of the first Lyapunov coefficient 7 (see 

(BSD). 
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